home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / coupling.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  12.7 KB  |  413 lines
  1. /* specfunc/coupling.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_errno.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_coupling.h>
  28.  
  29. #include "error.h"
  30.  
  31. inline
  32. static
  33. int locMax3(const int a, const int b, const int c)
  34. {
  35.   int d = GSL_MAX(a, b);
  36.   return GSL_MAX(d, c);
  37. }
  38.  
  39. inline
  40. static
  41. int locMin3(const int a, const int b, const int c)
  42. {
  43.   int d = GSL_MIN(a, b);
  44.   return GSL_MIN(d, c);
  45. }
  46.  
  47. inline
  48. static
  49. int locMin5(const int a, const int b, const int c, const int d, const int e)
  50. {
  51.   int f = GSL_MIN(a, b);
  52.   int g = GSL_MIN(c, d);
  53.   int h = GSL_MIN(f, g);
  54.   return GSL_MIN(e, h);
  55. }
  56.  
  57.  
  58. /* See: [Thompson, Atlas for Computing Mathematical Functions] */
  59.  
  60. static
  61. int
  62. delta(int ta, int tb, int tc, gsl_sf_result * d)
  63. {
  64.   gsl_sf_result f1, f2, f3, f4;
  65.   int status = 0;
  66.   status += gsl_sf_fact_e((ta + tb - tc)/2, &f1);
  67.   status += gsl_sf_fact_e((ta + tc - tb)/2, &f2);
  68.   status += gsl_sf_fact_e((tb + tc - ta)/2, &f3);
  69.   status += gsl_sf_fact_e((ta + tb + tc)/2 + 1, &f4);
  70.   if(status != 0) {
  71.     OVERFLOW_ERROR(d);
  72.   }
  73.   d->val = f1.val * f2.val * f3.val / f4.val;
  74.   d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val);
  75.   return GSL_SUCCESS;
  76. }
  77.  
  78.  
  79. static
  80. int
  81. triangle_selection_fails(int two_ja, int two_jb, int two_jc)
  82. {
  83.   return ((two_jb < abs(two_ja - two_jc)) || (two_jb > two_ja + two_jc));
  84. }
  85.  
  86.  
  87. static
  88. int
  89. m_selection_fails(int two_ja, int two_jb, int two_jc,
  90.                   int two_ma, int two_mb, int two_mc)
  91. {
  92.   return (   abs(two_ma) > two_ja 
  93.           || abs(two_mb) > two_jb
  94.       || abs(two_mc) > two_jc
  95.       || GSL_IS_ODD(two_ja + two_ma)
  96.       || GSL_IS_ODD(two_jb + two_mb)
  97.       || GSL_IS_ODD(two_jc + two_mc)
  98.           || (two_ma + two_mb + two_mc) != 0
  99.       );
  100. }
  101.  
  102.  
  103. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  104.  
  105. int
  106. gsl_sf_coupling_3j_e(int two_ja, int two_jb, int two_jc,
  107.                         int two_ma, int two_mb, int two_mc,
  108.             gsl_sf_result * result)
  109. {
  110.   /* CHECK_POINTER(result) */
  111.  
  112.   if(two_ja < 0 || two_jb < 0 || two_jc < 0) {
  113.     DOMAIN_ERROR(result);
  114.   }
  115.   else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
  116.           || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)
  117.      ) {
  118.     result->val = 0.0;
  119.     result->err = 0.0;
  120.     return GSL_SUCCESS;
  121.   }
  122.   else {
  123.     gsl_sf_result n1_a, n1_b, n3_a, n3_b;
  124.     gsl_sf_result d1_a, d1_b, d2_a, d2_b, d3_a, d3_b;
  125.     gsl_sf_result n1, n2, n3;
  126.     gsl_sf_result d1, d2, d3;
  127.     double norm;
  128.     double sign = (GSL_IS_ODD((two_ja - two_jb - two_mc)/2) ? -1.0 : 1.0);
  129.     int tk, tkmin, tkmax;
  130.     double sum_pos = 0.0;
  131.     double sum_neg = 0.0;
  132.     double phase;
  133.     int status = 0;
  134.     status += gsl_sf_fact_e((two_jc + two_ja - two_jb)/2, &n1_a);
  135.     status += gsl_sf_fact_e((two_jc - two_ja + two_jb)/2, &n1_b);
  136.     status += gsl_sf_fact_e((two_ja + two_jb - two_jc)/2, &n2);
  137.     status += gsl_sf_fact_e((two_jc - two_mc)/2, &n3_a);
  138.     status += gsl_sf_fact_e((two_jc + two_mc)/2, &n3_b);
  139.     status += gsl_sf_fact_e((two_ja + two_jb + two_jc)/2 + 1, &d1);
  140.     status += gsl_sf_fact_e((two_ja - two_ma)/2, &d2_a);
  141.     status += gsl_sf_fact_e((two_ja + two_ma)/2, &d2_b);
  142.     status += gsl_sf_fact_e((two_jb - two_mb)/2, &d3_a);
  143.     status += gsl_sf_fact_e((two_jb + two_mb)/2, &d3_b);
  144.  
  145.     if(status != 0) {
  146.       OVERFLOW_ERROR(result);
  147.     }
  148.  
  149.     n1.val = n1_a.val * n1_b.val;
  150.     n3.val = n3_a.val * n3_b.val;
  151.     d2.val = d2_a.val * d2_b.val;
  152.     d3.val = d3_a.val * d3_b.val;
  153.  
  154.     norm = sign * sqrt(n1.val*n2.val*n3.val)/sqrt(d1.val*d2.val*d3.val);
  155.  
  156.     tkmin = GSL_MAX(0, two_jb - two_ja - two_mc);
  157.     tkmax = GSL_MIN(two_jc - two_ja + two_jb, two_jc - two_mc);
  158.     
  159.     phase = GSL_IS_ODD((tkmin + two_jb + two_mb)/2) ? -1.0 : 1.0;
  160.  
  161.     for(tk=tkmin; tk<=tkmax; tk += 2) {
  162.       double term;
  163.  
  164.       status = 0;
  165.       status += gsl_sf_fact_e((two_jb + two_jc + two_ma - tk)/2, &n1);
  166.       status += gsl_sf_fact_e((two_ja - two_ma + tk)/2, &n2);
  167.       status += gsl_sf_fact_e(tk/2, &d1_a);
  168.       status += gsl_sf_fact_e((two_jc - two_ja + two_jb - tk)/2, &d1_b);
  169.       status += gsl_sf_fact_e((two_jc - two_mc - tk)/2, &d2);
  170.       status += gsl_sf_fact_e((two_ja - two_jb + two_mc + tk)/2, &d3);
  171.  
  172.       if(status != 0) {
  173.         OVERFLOW_ERROR(result);
  174.       }
  175.  
  176.       d1.val = d1_a.val * d1_b.val;
  177.  
  178.       term = phase * n1.val * n2.val / (d1.val * d2.val * d3.val);
  179.       phase = -phase;
  180.  
  181.       if(norm*term >= 0.0) {
  182.         sum_pos += norm*term;
  183.       }
  184.       else {
  185.         sum_neg -= norm*term;
  186.       }
  187.     }
  188.  
  189.     result->val  = sum_pos - sum_neg;
  190.     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  191.     result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin) * fabs(result->val);
  192.  
  193.     return GSL_SUCCESS;
  194.   }
  195. }
  196.  
  197.  
  198. int
  199. gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc,
  200.                         int two_jd, int two_je, int two_jf,
  201.             gsl_sf_result * result)
  202. {
  203.   /* CHECK_POINTER(result) */
  204.  
  205.   if(   two_ja < 0 || two_jb < 0 || two_jc < 0
  206.      || two_jd < 0 || two_je < 0 || two_je < 0
  207.      ) {
  208.     DOMAIN_ERROR(result);
  209.   }
  210.   else if(   triangle_selection_fails(two_ja, two_jb, two_je)
  211.           || triangle_selection_fails(two_ja, two_jc, two_jf)
  212.           || triangle_selection_fails(two_jb, two_jd, two_jf)
  213.           || triangle_selection_fails(two_jc, two_jd, two_je)
  214.      ) {
  215.     result->val = 0.0;
  216.     result->err = 0.0;
  217.     return GSL_SUCCESS;
  218.   }
  219.   else {
  220.     gsl_sf_result n1;
  221.     gsl_sf_result d1, d2, d3, d4, d5, d6;
  222.     double norm;
  223.     int tk, tkmin, tkmax;
  224.     double phase;
  225.     double sum_pos = 0.0;
  226.     double sum_neg = 0.0;
  227.     double sumsq_err = 0.0;
  228.     int status = 0;
  229.     status += delta(two_ja, two_jb, two_je, &d1);
  230.     status += delta(two_ja, two_jc, two_jf, &d2);
  231.     status += delta(two_jb, two_jd, two_jf, &d3);
  232.     status += delta(two_jc, two_jd, two_je, &d4);
  233.     if(status != GSL_SUCCESS) {
  234.       OVERFLOW_ERROR(result);
  235.     }
  236.     norm = sqrt(d1.val) * sqrt(d2.val) * sqrt(d3.val) * sqrt(d4.val);
  237.     
  238.     tkmin = locMax3(0,
  239.                    two_ja + two_jd - two_je - two_jf,
  240.                    two_jb + two_jc - two_je - two_jf);
  241.  
  242.     tkmax = locMin5(two_ja + two_jb + two_jc + two_jd + 2,
  243.                     two_ja + two_jb - two_je,
  244.             two_jc + two_jd - two_je,
  245.             two_ja + two_jc - two_jf,
  246.             two_jb + two_jd - two_jf);
  247.  
  248.     phase = GSL_IS_ODD((two_ja + two_jb + two_jc + two_jd + tkmin)/2)
  249.             ? -1.0
  250.         :  1.0;
  251.  
  252.     for(tk=tkmin; tk<=tkmax; tk += 2) {
  253.       double term;
  254.       double term_err;
  255.       gsl_sf_result den_1, den_2;
  256.       gsl_sf_result d1_a, d1_b;
  257.       status = 0;
  258.       
  259.       status += gsl_sf_fact_e((two_ja + two_jb + two_jc + two_jd - tk)/2 + 1, &n1);
  260.       status += gsl_sf_fact_e(tk/2, &d1_a);
  261.       status += gsl_sf_fact_e((two_je + two_jf - two_ja - two_jd + tk)/2, &d1_b);
  262.       status += gsl_sf_fact_e((two_je + two_jf - two_jb - two_jc + tk)/2, &d2);
  263.       status += gsl_sf_fact_e((two_ja + two_jb - two_je - tk)/2, &d3);
  264.       status += gsl_sf_fact_e((two_jc + two_jd - two_je - tk)/2, &d4);
  265.       status += gsl_sf_fact_e((two_ja + two_jc - two_jf - tk)/2, &d5);
  266.       status += gsl_sf_fact_e((two_jb + two_jd - two_jf - tk)/2, &d6);
  267.       
  268.       if(status != GSL_SUCCESS) {
  269.         OVERFLOW_ERROR(result);
  270.       }
  271.  
  272.       d1.val = d1_a.val * d1_b.val;
  273.       d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err;
  274.  
  275.       den_1.val  = d1.val*d2.val*d3.val;
  276.       den_1.err  = d1.err * fabs(d2.val*d3.val);
  277.       den_1.err += d2.err * fabs(d1.val*d3.val);
  278.       den_1.err += d3.err * fabs(d1.val*d2.val);
  279.  
  280.       den_2.val  = d4.val*d5.val*d6.val;
  281.       den_2.err  = d4.err * fabs(d5.val*d6.val);
  282.       den_2.err += d5.err * fabs(d4.val*d6.val);
  283.       den_2.err += d6.err * fabs(d4.val*d5.val);
  284.  
  285.       term  = phase * n1.val / den_1.val / den_2.val;
  286.       phase = -phase;
  287.       term_err  = n1.err / fabs(den_1.val) / fabs(den_2.val);
  288.       term_err += fabs(term / den_1.val) * den_1.err;
  289.       term_err += fabs(term / den_2.val) * den_2.err;
  290.  
  291.       if(term >= 0.0) {
  292.         sum_pos += norm*term;
  293.       }
  294.       else {
  295.         sum_neg -= norm*term;
  296.       }
  297.  
  298.       sumsq_err += norm*norm * term_err*term_err;
  299.     }
  300.  
  301.     result->val  = sum_pos - sum_neg;
  302.     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  303.     result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
  304.     result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val);
  305.  
  306.     return GSL_SUCCESS;
  307.   }
  308. }
  309.  
  310.  
  311. int
  312. gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc,
  313.                         int two_jd, int two_je, int two_jf,
  314.             int two_jg, int two_jh, int two_ji,
  315.             gsl_sf_result * result)
  316. {
  317.   /* CHECK_POINTER(result) */
  318.  
  319.   if(   two_ja < 0 || two_jb < 0 || two_jc < 0
  320.      || two_jd < 0 || two_je < 0 || two_jf < 0
  321.      || two_jg < 0 || two_jh < 0 || two_ji < 0
  322.      ) {
  323.     DOMAIN_ERROR(result);
  324.   }
  325.   else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
  326.           || triangle_selection_fails(two_jd, two_je, two_jf)
  327.           || triangle_selection_fails(two_jg, two_jh, two_ji)
  328.           || triangle_selection_fails(two_ja, two_jd, two_jg)
  329.           || triangle_selection_fails(two_jb, two_je, two_jh)
  330.           || triangle_selection_fails(two_jc, two_jf, two_ji)
  331.      ) {
  332.     result->val = 0.0;
  333.     result->err = 0.0;
  334.     return GSL_SUCCESS;
  335.   }
  336.   else {
  337.     int tk;
  338.     int tkmin = locMax3(abs(two_ja-two_ji), abs(two_jh-two_jd), abs(two_jb-two_jf));
  339.     int tkmax = locMin3(two_ja + two_ji, two_jh + two_jd, two_jb + two_jf);
  340.     double sum_pos = 0.0;
  341.     double sum_neg = 0.0;
  342.     double sumsq_err = 0.0;
  343.     double phase;
  344.     for(tk=tkmin; tk<=tkmax; tk += 2) {
  345.       gsl_sf_result s1, s2, s3;
  346.       double term;
  347.       double term_err;
  348.       int status = 0;
  349.       status += gsl_sf_coupling_6j_e(two_ja, two_ji, two_jd,  two_jh, tk, two_jg,  &s1);
  350.       status += gsl_sf_coupling_6j_e(two_jb, two_jf, two_jh,  two_jd, tk, two_je,  &s2);
  351.       status += gsl_sf_coupling_6j_e(two_ja, two_ji, two_jb,  two_jf, tk, two_jc,  &s3);
  352.       if(status != GSL_SUCCESS) {
  353.         OVERFLOW_ERROR(result);
  354.       }
  355.       term = s1.val * s2.val * s3.val;
  356.       term_err  = s1.err * fabs(s2.val*s3.val);
  357.       term_err += s2.err * fabs(s1.val*s3.val);
  358.       term_err += s3.err * fabs(s1.val*s2.val);
  359.  
  360.       if(term >= 0.0) {
  361.         sum_pos += (tk + 1) * term;
  362.       }
  363.       else {
  364.         sum_neg -= (tk + 1) * term;
  365.       }
  366.  
  367.       sumsq_err += ((tk+1) * term_err) * ((tk+1) * term_err);
  368.     }
  369.  
  370.     phase = GSL_IS_ODD(tkmin) ? -1.0 : 1.0;
  371.  
  372.     result->val  = phase * (sum_pos - sum_neg);
  373.     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  374.     result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
  375.     result->err += 2.0 * GSL_DBL_EPSILON * (tkmax-tkmin + 2.0) * fabs(result->val);
  376.  
  377.     return GSL_SUCCESS;
  378.   }
  379. }
  380.  
  381.  
  382. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  383.  
  384. #include "eval.h"
  385.  
  386. double gsl_sf_coupling_3j(int two_ja, int two_jb, int two_jc,
  387.                           int two_ma, int two_mb, int two_mc)
  388. {
  389.   EVAL_RESULT(gsl_sf_coupling_3j_e(two_ja, two_jb, two_jc,
  390.                                    two_ma, two_mb, two_mc,
  391.                                    &result));
  392. }
  393.  
  394.  
  395. double gsl_sf_coupling_6j(int two_ja, int two_jb, int two_jc,
  396.                           int two_jd, int two_je, int two_jf)
  397. {
  398.   EVAL_RESULT(gsl_sf_coupling_6j_e(two_ja, two_jb, two_jc,
  399.                                    two_jd, two_je, two_jf,
  400.                                    &result));
  401. }
  402.  
  403.  
  404. double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc,
  405.                           int two_jd, int two_je, int two_jf,
  406.                           int two_jg, int two_jh, int two_ji)
  407. {
  408.   EVAL_RESULT(gsl_sf_coupling_9j_e(two_ja, two_jb, two_jc,
  409.                                    two_jd, two_je, two_jf,
  410.                                    two_jg, two_jh, two_ji,
  411.                                    &result));
  412. }
  413.